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Abstract —In this paper, we propose an optimization frame¬ 
work to design a network of positive linear systems whose 
structure switches according to a Markov process. The optimiza¬ 
tion framework herein proposed allows the network designer to 
optimize the coupling elements of a directed network, as well as 
the dynamics of the nodes in order to maximize the stabilization 
rate of the network and/or the disturbance rejection against an 
exogenous input. The cost of implementing a particular network 
is modeled using posynomial cost functions, which allow for a 
wide variety of modeling options. In this context, we show that 
the cost-optimal network design can be efficiently found using 
geometric programming in polynomial time. We illustrate our 
results with a practical problem in network epidemiology, namely, 
the cost-optimal stabilization of the spread of a disease over a 
time-varying contact network. 


I. Introduction 


T he intricate structure of many biological, social, and 
economic networks emerges as the result of local inter¬ 
actions between agents aiming to optimize their utilities. The 
emerging networked system must satisfy both structural and 
functional requirements, even in the presence of time-varying 
interactions. An important set of functional requirements is 
concerned with the behavior of dynamic processes taking place 
in the network. For example, most biological networks emerge 
as the result of an evolutive process that forces the network 
to be stable and robust to external perturbations. 

Among networks of dynamical systems, those consisting 
of positive systems (i.e., the state variables are nonnegative 
quantities provided the initial state and the inputs are non¬ 
negative are of particular importance. Positive systems 
arise naturally while modeling systems in which the variables 
of interest are inherently nonnegative, such as concentration 
of chemical species [12;|, information rates in communication 
networks (21) , sizes of infected populations in epidemiol¬ 
ogy fl7) , (19[ , and many other compartmental models |3). 
Due to its practical relevance, many control-theoretical tools 
have been adapted to the particular case of positive linear 
systems, such as the characterization of stability via diagonal 
Lyapunov functions (2j, the bounded real lemma |22|, and 
integral linear constraints for robust stability (7), to mention 
a few (we point the reader to for a thorough exposition 
on positive systems). 

Most of the methods mentioned above, however, focus on 
state-feedback control of positive linear systems. In contrast. 


there are many situations in which state-feedback is not a 
feasible option. A particular example is the stabilization of 
a viral spreading process in a complex contact network ED’ 
& where it is not feasible to obtain reliable measurements 
about the state of the spread in (almost) real time. In this 
and similar situations, the control problem is better posed as 
the problem of tuning the dynamics of the nodes and the 
coupling elements of the network in order to guarantee a 
stable dynamics. Furthermore, existing control methods mostly 
focus on the analysis of time-invariant systems; thus, they 
do not provide effective tools for designing time-switching 
topologies. 

The aim of this paper is to propose a tractable optimization 
framework to design networked positive linear systems in the 
presence of time-switching topologies. We model the time 
variation of the network structure using Markov processes, 
where the modes of the Markov process correspond to dif¬ 
ferent network structures and the transition rates indicate the 
probability of switching between topologies. In this context, 
we consider the problem of designing both the dynamics of 
the nodes in the network, as well as the structure of the 
elements coupling them in order to optimize the dynamic 
performance of the switching network. We model the cost of 
implementing a particular network dynamics using posynomial 
cost functions, which allow for a wide range of modeling op¬ 
tion [,5j. We then propose an efficient optimization framework, 
based on geometric programming [5), to find the cost-optimal 
network design that maximizes the stabilization rate and/or the 
disturbance rejection against exogenous signals. To achieve 
this objective, we develop new theoretical characterizations of 
the stabilization rate and the disturbance attenuation of positive 
Markov jump linear systems which are specially amenable in 
the context of geometric programming. 

The paper is organized as follows. In Section [II] we in¬ 
troduce elements of graph theory, Markov jump linear sys¬ 
tems, and geometric programming used in our derivations. In 
Section III we model the dynamics of randomly switching 
network using Markov jump linear systems and rigorously 
state the design problems under consideration. Then, we 
propose geometric programs to efficiently find the cost-optimal 
network design for stabilization and disturbance attenuation, in 
Sections |IV] and [V] respectively. We illustrate our results with 
a relevant epidemiological problem, namely, the stabilization 


of a viral spread in a switching contact network, in Section VI 
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II. Mathematical Preliminaries 
In this section, we introduce the notation and some re¬ 
sults needed in our derivations. We denote by R. and R + 





2 


the sets of real and nonnegative numbers, respectively. The 
set {1, ...,7V} is sometimes denoted by [TV]. We denote 
vectors using boldface letters and matrices using capital letters. 
When x £ R" is nonnegative (positive) entrywise, we write 
x > 0 (x > 0, respectively). The 1-norm on R” is denoted by 
Hi = £r=iM- The n-vector and n x m matrix whose en¬ 
tries are all one are denoted by l n and l nxm , respectively. We 
denote the identity matrix by I. The (z, j)-entry of a matrix A 
is denoted by A t j. The nxm matrix whose (z, j)-entry equals 
dij is denoted by [aij]j e [n] je[m]- or simply by when n 

and m are clear from the context. A square matrix is Hurwitz 
stable if the real parts of its eigenvalues are all negative. A 
square matrix is Metzler if its off-diagonal entries are nonneg¬ 
ative. The Kronecker product | 6 j of the matrices A and B is 
denoted by A 0 B. Following fil, we define the generalized 
Kronecker product of an n\ x ri 2 matrix A = [dij\i.-j and the 
set of mi x m 2 matrices }ie[ni]je[n 2 ] as the following 

(rtimi) x (ri 2 m 2 ) block matrix 


A 0 {Bij}ij 


aiiBn 

(t-n \ 1 B ni 1 


Ctlri2 C^ln2 


tnin2-Unin2_ 


The direct sum of the square matrices A\, ..., A n , denoted by 
©'Li A t , is defined as the block-diagonal matrix containing 
the matrices A\,..., A n as its diagonal blocks. 


A. Basic Graph Theory 

A weighted, directed graph (also called digraph) is defined 
as the triad Q = (V, 8, W), where V = (iq,..., v n } is a set 
of n nodes, £ C V x V is a set of ordered pairs of nodes called 
directed edges, and the function W : £ — t ( 0 ,oo) assigns 
positive real weights to the edges in £. By convention, we 
say that ( Vi,Vj) is an edge from Vj pointing towards Vj. The 
adjacency matrix of a weighted and directed graph Q is an 
nxn matrix [a,;,],;., defined entry-wise as = W((vi,Vj )) if 
edge ( Vi,Vj ) £ £, and aij = 0 otherwise. Thus, the adjacency 
matrix of a graph is always nonnegative. Conversely, given a 
nonnegative matrix A, we can associate to it a directed graph 
whose adjacency matrix is A. 


B. Positive Markov Jump Linear Systems 
Consider the linear time-invariant system 

x = Ax + Bw, z = Cx + Dw, (1) 

where x(f) £ R", w(t) £ R s , z(t) S R r , and A, B, C, and D 
are real matrices of appropriate dimensions. The system 0 is 
usually denoted by the quadruple ( A , B , G, D). We say that 
the system 0 is positive if x(f) > 0 and z(f) > 0 for all 
t > 0 provided that x( 0 ) = Xo > 0 and w(f) > 0 for every 
t > 0. It is well known m that the system 0 is positive if 
and only if A is Metzler and B, C, and D are nonnegative. 

Let a = {er(t)} t >o be a time-homogeneous Markov process 
with the state space {1,... ,M}. Let fl = [7R ? L je[Af] denote 
the infinitesimal generator of a, i.e., assume that the transition 
probability of o is given by Pr(cr(i + h) = j \ a(t) = i) = 
ttijh + o{h) (i ^ j) and Pr(cr(t + h) = i | aft) = i) = 


l+Ttiih + o(h) for all z, j £ [TIL] and t, h > 0. We will assume 
that the Markov process is irreducible. A Markov jump linear 
system | 8 ] is described as the following stochastic differential 
equations: 

s . f *(*) = A a(t)*(t) + B a{t) w{t), 

' \z(f) = C CT(t) x(f) + D a{t) w(t). 

We assume that the initial states x(0) = x 0 and rr(0) = do are 
constants. We say that the Markov jump linear system E is 
positive if the linear time-invariant systems (Ai, Bi,Ci, Di) 
are positive for all i £ [M], We say that E is internally 
mean stable (or simply mean stable ) if there exist a > 0 
and A > 0 such that, if w(f) = 0 for every t > 0, then 
£ , [||x(t)|| 1 ] < ae _At ||x 0 ||i for all x 0 , cr 0 , and t > 0 , where 
E denotes the expectation operator. The supremum of A such 
that this inequality holds for all Xo, cr o, and t is called the 
exponential decay rate of a mean stable E. The following 
proposition describes necessary and sufficient conditions for 
mean-stability of positive Markov jump linear systems. 

Proposition 1 (7[?| J, |75]/): If the Markov jump linear sys¬ 
tem E is positive, then the following conditions are equivalent: 

1) E is mean stable. 

2) The matrix A = II T 0 I + ©]L | Ai is Hurwitz stable. 

3) There exist positive vectors Vi,..., vm G R" such that 

A i + £j=i 7r R' v J < 0 for every * € [M], 

We let /Tx (M.+, ) denote the space of Lebesgue in- 
tegrable functions on R + taking values in R!}- F° r an Y 
/ e £i(R+,R"), we define ||/|| £l = / 0 °°||/(t )|| 1 dt. The 
definition below extends the concept of Ei-stability (given in, 
e.g., 0 ) to Markov jump linear systems: 

Definition 1: We say that E is C\-stable if there exists 7 > 
0 such that, for all cro and w £ £i(R + ,RL), we have that 
E[z] £ £i(R + ,RIj_) and ||.E[z]||, Cl < 7 ||w ||£ 1 when x 0 = 0. 
If E is Ex-stable then its C\-gain, denoted by [E|| 1; is defined 
by ||E||x =sup we£i(R+iR »_ ) (||L;[z]|| z:i /||w|| £l ). 

For the particular case when E is time-invariant, the next 
characterization of the C \ -gain is available: 

Proposition 2 mm : Assume that the linear time- 
invariant system 0 is positive and let 7 > 0 be arbitrary. 
The following statements are equivalent: 

1) The system 0 is stable and its Ei-gain is less than 7 . 

2) A is Hurwitz stable and {D — CA~ 1 B) — 7 IJ" < 0. 

3) There exists a positive vector v £ R n satisfying inequal¬ 
ities v T A + t^C < 0 and v T B + i^D < 7 l}L 


C. Geometric Programming 

The design framework proposed in this paper depends on a 
class of optimization problems called geometric programs (5j. 
Let x-i , ..., x n denote n real positive variables and define the 
vector variable x = (27, ... ,x n ). We say that a real-valued 
function g(x) is a monomial function ( monomial for short) 
if there exist c > 0 and ai,... ,a n £ R such that <?(x) = 
cx • • • x an . Also we say that a real-valued function /(x) is 
a posynomial function (posynomial for short) if it is a sum of 
monomial functions of x. For information about the modeling 
power of posynomials, we point the reader to |5j. Given 
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positive integers p, q and a collection of posynomials /o(x), 
..., / p (x) and monomials <n(x), .... g q {x), the optimization 
problem 

minimize /o(x) 

X 

subject to fi (x) < 1 , i = l,...,p, (3) 

5 j(x) = 1 , j = l,...,q, 

is called a geometric program (GP) in standard form. Although 
geometric programs are not convex, they can be efficiently 
converted into a convex optimization problem and efficiently 
solved using, for example, interior-point methods (see ©, for 
more details on GP). We quote the following result from tm 
Section 10.4] on the computational complexity of solving GP: 

Proposition 3: The geometric program Q can be solved 
with computational cost 0((k+m) 1 ^ 2 (mk 2 +k 3 +n 3 )), where 
m = p + 2q and k is the maximum of the numbers of mono¬ 
mials contained in each of posynomials /o(x), ..., / p (x). 

Geometric programs can also be written in terms of matrices 
and vectors, as follows. Let X-[,..., X n be matrix-valued 
positive variables and define X = {X -\,..., X n ). We say that a 
real matrix-valued function F(X) is a monomial (posynomial) 
if each entry of F(X) is a monomial (respectively, posyno¬ 
mial) in the entries of the matrix-valued variables Xi, ..., 
X n . Then, given a scalar-valued posynomial /o(X), matrix- 
valued posynomials Fi(X), ..., F p (X), and matrix-valued 
monomials Gi(X), ..., G q (X), we call the optimization 
problem 

minimize /o(X) 

subject to Fj(X) <1, i = 1,... ,p, 

Gj(X) = 1, j = 1,... ,q, 

a matrix geometric program, where the right-hand side of the 
constraints are all-ones matrices of appropriate dimensions. 
Notice that one can easily reduce a matrix geometric program 
to a standard geometric program by dealing with the matrix¬ 
valued constraints entry-wise. 

The next lemma will be used in Sections [IV] and [V] to prove 
our main results. 

Lemma 1: Let X fh X \, ..., X n be matrix-valued pos¬ 
itive variables. Assume that Xq is RP X 9-valued, and let 
F(X) be an R px<? -valued posynomial in the variable X = 
(Xq, X-[,, X n ). Then, there exists a matrix-valued posyno¬ 
mial F(X) such that F(X) < Xq if and only if F(X) < t pxq 
for every X. 

Proof: Define the function F of X component-wise as 
F(X) k/ = (X 0 ) k jF(X) kt t for all k G \p\ and l G [q\. Then, 
F(X) is a matrix-valued posynomial. Then, since X [} > 0, the 
constraint F(X) < X 0 is equivalent to F(X) < t pxq . ■ 


III. Switched Network of Positive Linear Systems 

In this section, we describe the dynamics of the network of 
positive linear systems under consideration and state relevant 
design problems. 


A. Network Dynamical Model 

Consider a collection of N linear time-invariant subsystems 


= F k x k + G k u fc , y k = -fffcXfc, 


(4) 


where k G V = {1 x fc (f) G y k (t) G 

and u k (t) G R. qk . Assume that these subsystems are lin¬ 
early coupled through the edges of a time-varying, weighted, 
and directed graph having a time-dependent adjacency ma¬ 
trix I\(t) = {K k i{t)\ k ,i^v■ The coupling between subsystems 
is modeled by the following set of inputs: 


u k(t) = ^K M (t)T (k/) y t {f), 


(5) 


£ev 


where k G V and I, 


(M) 


oXpf 


is the so-called inner 


coupling matrix, which indicates how the output of the f-th 
subsystem influences the input of the fc-th subsystem. The 
inner-coupling matrix can be interpreted as a matrix-valued 
weight associated to each directed edge of the graph. 

In the rest of the paper, we consider dynamic networks 
satisfying the following positivity assumption: 

Assumption 1: System 0 is positive and the matrix T(M) 
is nonnegative for all k, £ G V. 

The above assumption is satisfied for many networked 
dynamics where the state variables are nonnegative. For ex¬ 
ample, the dynamics of many models of disease spreading 
in networks | [23] can be described as a coupled network of 
positive systems, where the state variables represent prob¬ 
abilities of infection. We will present the details of this 


model in Section VI Other examples of practical relevance 
are chemical networks m transportation networks ( 20 ), and 
compartmental models 0 . 

In this paper, we consider networks of positive linear 
systems in which the network structure switches according 
to a Markov process, as indicated below: 

Assumption 2: There exist N x N nonnegative matri¬ 
ces K\, ..., Km and a time-homogeneous Markov pro¬ 
cess a = {cr(f)}t>o taking its values in [M\ such that 
K(t) = K a{t) for every t > 0. 

Under the above assumptions, the dynamics of the net¬ 
work of subsystems in 0 coupled through the law 0 
forms a positive Markov jump linear system, as described 
below. Let us define the ‘stacked’ vectors x = (x^~,..., x|.) 1 , 
u = (u^,.. • ,u^) T , and y = (y^, ■ ■ •,y^) T . Then, the dy- 
namics in 0 can be written by the differential equations x = 
(©f=iTfe)x+ (®f = 1 G fe )u and y = (Q) k=1 H k )x. Simi¬ 
larly, (5]i can be written using the generalized Kronecker prod¬ 
uct (defined in Section [U i as u(t) = (if(f)®{r (M } k ,e)y(t). 
Thus, the global network dynamics, denoted by Af, can be 
written as the following positive Markov jump linear system: 

TV : x(i) = A aW x(t), 

where the matrices Ai (i = 1,..., M) are given by 


N 


N 


* = ® F k + © Gfcj f Ki<8>{T(k,e)}k,t 0 H k ). ( 6 ) 

k =l l ' ^ ^ l 


N 
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B. Network Design: Cost and Constraints 

In this paper, we propose a novel methodology to simulta¬ 
neously design the dynamics of each subsystem (characterized 
by the matrices F k , G k , and II k for all k £ V) and the 
weights of the edges coupling them (characterized by the 
inner-coupling matrices r (M ) for k, t £ V). Our design 
objective is to achieve a prescribed performance criterion 
(described below) while satisfying certain cost requirements. 
In our problem formulation, the time-variant graph structure 
changes according to a Markov process er that we assume 
to be an exogenous signal. In other words, we assume the 
Markov process ruling the network switching is out of our 
control. We refer the readers to GD where, for a model of 
spreading processes, the authors study an optimal network 
design problem over networks whose structure endogenously 
changes by reacting to the dynamic state of agents. On the 
other hand, we assume that we can modify the dynamics of 
the subsystems and the matrix weights of each edge to achieve 
our design objective. 

We specifically assume that, for all k and t in V, there is 
a cost function associated to its coupling matrix r {k,€>- We 
denote this real-valued edge-cost function by 4>(k,i)- This cost 
function can represent, for example, the cost of building an 
interconnection between two subsystems. We remark that, in 
fact, we do not need to design T(M) for all k and £ by the 
following reason. For each i £ [M\, let Qi = 
denote the weighted directed graph having the adjacency 
matrix K ,, and consider the union £ = U ie[M]^i °f the 
possible directed edges. Then we see that, if (k, £) ^ £, then 
the Markov jump linear system A F under our consideration is 
independent of the value of T \ k ,t)- Therefore, throughout the 
paper, we focus on designing the coupling matrices T e only for 
e £ £. Similarly, our framework allows us to associate a cost to 
each subsystem in the network. For each k £ V, we denote the 
cost of implementing the k -th subsystem by fk{F k ,G k , H k ). 
Thus, the total cost of realizing a particular network dynamics 
is given by 

= + 0e(r e ). (7) 

fcGV e6f 

In practice, not all realizations of the subsystem 
(Fk,Gk, Hk) and the coupling matrix r e are feasible. We 
account for feasibility constraints on the design of the k -th 
subsystem via the following set of inequalities and equalities: 

9k, P (Fk,Gk,Hk) < 1, hX >q (F k ,G k ,H k )= 1, (8) 

where p £ [s*;], q £ [ffc], and g kp and h k are real functions 
for each k £ V. Similarly, for each edge e £ £, we account 
for design constraints on the coupling matrix F e via the 
restrictions: 

9e, P (T e ) < 1, /4(r e ) = l, (9) 

where p £ [s e ], q £ [t e ], and g^ and hf q are real functions 
for each e £ £. 


TABLE I: Design problems under consideration. 



Performance-Constr. 

Budget-Cons tr. 

Stabilization Rate 

Problem I-A 

Problem I-B 

Disturbance Atten. 

Problem II-A 

Problem II-B 


C. Network Design: Problem Statements 

We are now in conditions to formulate the network design 
problems under consideration. Our design problems can be 
classified according to two different criteria. The first criterion 
is concerned with the dynamic performance. In this paper, 
we limit our attention to two performance indexes: (I) stabi¬ 
lization rate (considered in Section |Tv)), and (II) disturbance 
attenuation (considered in Section 0 . The second criterion is 
concerned with the design cost. According to this criterion, 
we have two types of design problems: (A) performance- 
constrained problems and ( B ) budget-constrained problems. 
In a budget-constrained problem, the designer is given a fixed 
budget and she has to find the network design to maximize 
a dynamic performance index (either the stabilization rate 
or the disturbance attenuation). In a performance-constrained 
problem, the designer is required to design a network that 
achieves a given performance index while minimizing the 
total cost of the design. Combinations of these two criteria 
described above result in four possible problem formulations, 
represented in Table |T| 

We describe Problems I-A and I-B in more rigorous terms 
in what follows (Problems II-A and II-B will be formulated 
in Section |V): 

Problem I-A (Performance-Constrained Stabilization): 
Given a desired decay rate A > 0, design the nodal dynamics 
{F k ,G k ,H k } keV and the coupling matrices {T e } e&£ such 
that the global network dynamics A f achieves an exponential 
decay rat(Q greater than A at a minimum implementation 
cost R, defined in 0 - while satisfying the feasibility 
constraints ([ 8 ]i and ([9]). 

Problem I-B (Budget-Constrained Stabilization): Given an 
available budget R > 0, design the nodal dynamics 
{F k ,G k ,H k } keV and the coupling matrices {r e } e6f in or¬ 
der to maximize the exponential decay rate A of Af while 
satisfying the budget constraint R < R, and the feasibility 
constraints <[ 8 ]> and ([9]). 

Remark 1: The optimal solutions of Problems I-A and I-B 
are inversely related, as explained below. Let R*( A) (A *(R)) 
be the optimal solution R of Problem I-A (A of Problem I-B, 
respectively). For simplicity in our presentation, let us assume 
that both R* and A* are strictly increasing functions and 
have nonempty domains. Also, suppose that the composi¬ 
tions A* o R* and R* o A* are well-defined. Then, from 
the definitions of A* and Rf, we have A*(i?*(A)) > A and 
R*(X*(R)) < R for all possible values of A and R. From 
these inequalities, we obtain A*(f?*(A*(F?))) > X*(R) and 
X*(R*(X*(R))) < A*(/?), respectively. We therefore have 
X*(R*(X*(R))) = X*(R). This yields that A* o R* equals 
the identity since A* is strictly increasing. In the same way, 

1 We remark that the exponential decay rate of N is well-defined, since N 
is a Markov jump linear system. 
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we can show that R* o A* is the identity. Hence, R* and A* 
are the inverses of each other. 

In the next section, we proceed to present an optimization 
framework to solve Problems I-A and I-B. In Section [V] we 
shall extend our results to Problems II-A and II-B, which 
we call the Performance-Constrained and Budget-Constrained 
Disturbance Attenuation problems, respectively. 


IV. Optimal Design for Network Stabilization 


The aim of this section is to present geometric programs 
in order to solve both the Performance- and the Budget- 
Constrained Stabilization problems, which are one of the main 
contributions of this paper. In what follows, we place the 
following assumption on the cost and constraint functions: 

Assumption 3: 

1) For each e £ £, the functions </> e and gf (p £ [s e ]) are 
posynomials, while hf q (q £ [t e ]) is a monomial. 

2) For each k £ V, there exists a real and diag¬ 
onal matrix £ W nkXnk such that, the func¬ 

tions f k (F k ,G kl H k ) = f k (F k - A k ,G k ,H k ) and 
g% p (F k ,G k ,H k ) = g k p (F k - A k ,G k ,H k ) (p £ [a*]) 
are posynomials, while the function (F k , G k , H k ) = 
h X, q (Fk - A k ,G k , H k ) ( q £ [f fe ]) is a monomial. 

Remark 2: Since the state matrix of a positive system is 
Metzler, it can contain negative diagonal elements. On the 
other hand, the decision variables of a geometric program 
must be positive. Therefore, negative diagonal entries cannot 
necessarily be directly used as decision variables in our 
optimization framework. The above assumption will be used 
to overcome this limitation and will allow us to design the 
diagonal elements by a suitable change of variables. 

In order to transform Problem I-A into a geometric program, 
we need to introduce the following definitions. Let 

5 = max(-7r ii ) + max (A fe )jj. (10) 

ie[M] kev,ie[n k ] 

Then, for every i E [M], define the nonnegative matrix 
Pi = (fu + S)I — ©^! A k- Also, for M nfcXnfe -valued pos¬ 
itive variables Fk (k E V), we define 


N 


k =1 


N 


k =1 


4 = ® F k + © G k J ( Ki g {r (M) } M ) ( (J) H k 


N 


k—1 


The next theorem shows how to efficiently solve the 
Performance-Constrained Stabilization problem via geometric 
programming: 

Theorem 1: The network design that solves Problem I-A 
is defined by the set of subsystems {(F£, G* kl H£)} k€ .v with 
Fj* = F 1 * _ Afc, and the coupling matrices {r*} e6 f, where 
the starred matrices are the solutions to the following matrix 


geometric program: 


minimize 

{Fk,Gk,Hk}k£V, 
{Te}eE£,{v;}i(=[M] 


Y, fk(F k ,G k ,H k ) + ^MFe) 

keV eGS 


(11a) 


subject to vj (A t + Pi + XI) + 7rij-vT < r)v[, (1 lb) 

g k , P {F k ,G k , H k ) < 1, (11c) 

hl q (F k ,G k ,H k ) = 1, (lid) 

©■ ( lle > 


Remark 3: The optimization program is, in fact, a 
matrix geometric program. The cost function in ( |1 1 a[ > is a 
posynomial under Assumption [3] Also, the set of design 
constraints < fTTc> (JTTeJ) are valid posynomial inequalities and 
monomial equalities. Furthermore, the constraint in ( |1 lb[ ) is a 
matrix-posynomial constraint by Lemma [I] Also the definition 
of S in ( [TO} ensures that the matrix Pi is nonnegative. 

Remark 4: Standard GP solvers cannot handle strict in¬ 
equalities, such as ( | 1 lb[ >. In practice, we can overcome this 
limitation by including an arbitrary small number to relax the 
strict inequality into a non-strict inequality. 

Before we present the proof of Theorem [T| we need to 
introduce the following corollary of Proposition |T| 

Corollary 1: Assume that the Markov jump linear system E 
defined in © is positive. Then, E is mean stable with an 
exponential decay rate greater than A > 0, if and only 
if there exist positive vectors vi,..., vm £ R n such that 

v 7 A i + 12jL i V J + < 0 for every * € [M], 

Proof: Let us assume that E is mean stable and its 
exponential decay rate is greater than A > 0. Then, the Markov 
jump linear system x(f) = (^(t) + A/)x(i) is mean stable. 
Therefore, by Proposition [T] we can find positive vectors 
Vj £ R" such that (Ai + XI) + < 0- This 

proves the necessity of the condition in the corollary. The 
sufficiency part can be proved in a similar way. ■ 

Let us prove Theorem [T] 

Proof of Theorem [7} By Corollary |T| the Performance- 
Constrained Stabilization problem is equivalent to the follow¬ 
ing optimization problem: 


minimize 'V] f k (F k , G k , H k ) + V] <j> e (F e ) (12a) 

{F k ,Gk,H k } kev . 

rp \ \ fcGV eGt 

M 

subject to vjA; + Y, TTjjvJ + AvJ < 0, Vj > 0, (12b) 

i=i 

© and 0. (12c) 

Notice that A,, defined in ©, can have negative diagonal 
entries, since F k is a Metzler matrix. If this is the case, 
the constraint in ( | 12b| > cannot be written as a posynomial 
inequality. To overcome this issue, for each k £ V, we use 
the following transformation 


F k = F k + A*;, (13) 


where A k is given in Assumption 0 In fact, if the 
triple ( F kl G k ,H k ) is a feasible solution of then 

f k (F k , G kl H k ) = f k (F k ,G kl H k ) is well-defined, which 
shows that F k is positive because f k is a posynomial. 
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Then, we show that ( | 12b| > is equivalent to ( |1 lb[ ). Noting 
that the transformation ( [T3 ] i yields A, = Ai + © fc=1 Afc, we 
can rewrite the constraint ( | 12 b| > as vjAi — wj 0 j V =1 At + 
ttuvj + + AVj 1 < 0 . Adding 6vJ to both 

sides of the above inequality, we obtain ( |1 lb| i. Also, the 
equivalence between the objective functions ( | 12 a| ) and (JTTaJ) 
is obvious from their definitions. In the same way, we can 
observe that the constraints and © are equivalent to 
the constraints ( |1 lc| > — ( | 11 e[ ). Therefore, we conclude that the 
optimization problems 0 and 0 are equivalent. Notice that 
the constraint v* > 0 is omitted in ( fTTj ) because ( fTT| is stated 
as a geometric program and, therefore, all of its variables are 
assumed to be positive. Finally, as mentioned in Remark [3] 
the optimization problem in 0 is indeed a matrix geometric 
program. ■ 

Theorem]!] allows us to find the cost-optimal network design 
to stabilize the system at a desired decay rate, assuming the 
feasible set defined by ( |1 lb| )-( fTTe} is not empty. Similarly, the 
following theorem introduces a geometric program to solve the 
Budget-Constrained Stabilization problem: 

Theorem 2: The network design that solves Problem I-B is 
described by the set of subsystems {(F£, G k , H k )} ke y with 
F£ = F£ — A^,, and the coupling matrices {r*} e6 £, where 
the starred matrices are the solutions to the following matrix 
geometric program: 

minimize 1 / A 

{r e }ee£:T v i}ie[M] 

subject to 0 ]>- 0 J, (14) 

fk(h, Gk,H k ) + 0e(r e ) < R. 

fcGV eGf 

Proof: By Corollary |T] we can formulate the Budget- 
Constrained Stabilization problem as the following optimiza¬ 
tion problem: 


maximize A (15) 

A ,{F k ,G k ,H k } k ev, 

(re}ee£,{ V i}ie[M] 

subject to <| 8 ). ©• (T2§, Vi > 0, A > 0, and R < R. 


Notice that maximizing A is equivalent to minimizing 1/A, 
since A > 0. The rest of the proof is similar to the proof 
of Theorem |T| and we limit ourselves to remark the main 
steps. First, we use the transformation in | [T3| to show that 
the optimization problem 0 is equivalent to 0 . We can 
show that the optimization problem ( fl4| is a matrix geometric 
program by following the explanation in Remark [3] ■ 

Remark 5: Using a discrete-time counterpart of Proposi¬ 
tion [7] one could easily establish discrete-time analogues 
of Theorems [I] and [2] in a straightforward manner (see, 
e.g., the stability characterizations of discrete-time positive 


Markov jump linear systems in [25 Theorem 1] and 115 
Theorem 3.4]). 


V. Optimal Design for Disturbance Attenuation 
In Section [IV] we have introduced an optimization frame¬ 
work, based on geometric programming, to solve both the 
Performance- and Budget-Constrained Stabilization problems 


described in Subsection |III-C| In this section, we extend 
our analysis to design networks from the point of view of 
disturbance attenuation. We consider the following collection 
of linear time-invariant subsystems: 

5 tfc — FfcXfc -p GfcpWfc T Gfc, 2 tifc, 

Zfc = #fc,lXfc + Jfc,uWfc + Jfc,12Ufc, (16) 

y k = -fffc, 2 x fc + 4,21 w fc, 

where k £ V. The signals w k (t) £ R Sfc and z k (t) £ are, 
respectively, the disturbance input and the performance output. 
We assume that subsystems are linearly interconnected through 
the law ([5]), and that the following positivity assumption holds: 

Assumption 4: For every k £ V, F k is Metzler and G k i, 
Gfc, 2 , H k> i, fTfc j2 , Jfc.n, Jk, 12 , and J fc)2 i are nonnegative. 

Using the generalized Kronecker product, we can rewrite 
the dynamics of the network Af of subsystems in ( fT~ 6 j i, coupled 
according to 0 . as a Markov jump linear system ([2]) with the 
following coefficient matrices: 


N 


N 


N 


A i ~ © F k + ( ® } k,z) 

k=l ^k—1 ' 'k =1 

N / iV v / N 

* = ® Gk,i + © Gfc, 2 1 (Ki®{r ik/) } kt t) © Jk, 


k =1 
N 


k =1 
N 


- © HkA + ( © 

k—1 ^k =1 ' i ' 

N , N v / N \ 

A = ® Jk,ll + © Jk, 12J {Ki®{F( k> e.)} k ,f) © Jk,2l \ • 

k=l ^fc=l ' 'k=l 


fc =1 
, N 


Therefore, the C \-gain (see Definition [TJ of the switched 
network of linearly coupled subsystems is well defined. We 
can then rigorously state Problem II-A, as follows: 

Problem II-A (Performance-Constrained Disturbance Attenuation): 
Given a desired C \ -gain, denoted by 7 > 0 , design 

the nodal dynamics {F k ,G kl H k } ke x> and the coupling 
matrices [ 1 / }«e£ such that the global network dynamics 
Af is mean stable and its C \-gain is less than 7, while the 
implementation cost R defined in ([7]i is minimized and the 
feasibility constraints 0 and © are satisfied. 

The next theorem, which gives the second main contribution 
of this paper, shows that this problem can be solved via 
geometric programming: 

Theorem 3: The network design that solves Problem II-A 
is defined by the set of subsystems {(F£, G k , Hf)} ks .\> with 
Ff = F£ — Afc, and the coupling matrices {r/} e6 £, where 
the starred matrices are the solutions to the following matrix 
geometric program: 


minimize 55 fk{F k ,G kl H k ) (r e ) 

\{Fk,G k ,Hk} kev , keV ee£ 

t i e>ee£)t v if<S[M] 

subject to vj(Ai + Pi) + 55 7 r ij v T + ^G t < SvJ, 
vjBi + tjDi < 7 IJ, 

(fTTc]-(fTTe|). 
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A. Proof 

The proof of Theorem [3] is based on the following theorem, 
which reduces the analysis of the C \ -gain of a general positive 
Markov jump linear system to a linear prograrr0 

Theorem 4: Assume that the Markov jump linear system E 
defined in 0 is positive. For any 7 > 0, E is internally 
mean stable and ||E||i < 7, if and only if there exist positive 
vectors vi,..., vju £ R n such that 


Ai 


M 

E 

i=i 




-j- l r Ci < 0, v i Bi + t r D, 


< 7ll, 


(17) 


for every i £ {1,..., M}. 

For the proof of Theorem |4j we introduce {0,1}-valued 
stochastic processes £ 1 , defined as £j(f) = 1 if 

<r(t) = i and £j(i) = 0 otherwise for every i and t, and define 
the vector £ = (£ lt ..., £m) T - Following the steps in the proof 


of 115 Proposition 5.3], one can prove that: 

+ B(E[(] <g> w), 


j t m ® *]= ae\z ( 

E[i ® z] = CE[£ ( 


D(E[£] ® w), 


(18) 


where A is defined in Proposition 0 B = ©fI T B u C = 
©^1 ©j and D = ©^ B, . We can then prove the following 
useful lemma: 

Lemma 2: If E is positive and internally mean stable, Xo = 
0 , and w(t) >0 for every t > 0 , then 



£[£ ® z \dt = (D 


CA~ l B) / (E[£]®w)dt. (19) 


Proof: Assume that E is internally mean stable and set 
xo = 0. Then, A is Hurwitz stable by Proposition |T| From 
( [T 8 ] >, it follows that £[£(t) ®z(t)] = f 0 Ce A ^~ s ^ .B(£[£(s)]® 
w(s)) ds + D(E[^(t)] ® w(t)). Since A is Hurwitz stable and, 
also, the functions Ce A ( t ~ s ' ) B, E[£ ® z], and E[£\ g) w are 
nonnegative, integrating this equation from 0 to 00 with respect 
to t completes the proof. ■ 

We also prove the following lemma, which can provide 
alternative expressions for HwH^ and ||£[z]||£ 1 : 

Lemma 3: Assume that E is positive. If Xo = 0 and w (t) > 
0 for every t > 0 , then 


pOO pOO 

/ || w (f)||i dt = tj ls (£[£] <g> w) eft, ( 20 ) 

Jo Jo 

p oo pOC 

/ \\ E Ht)]\\i dt = tjfr / E[f®z]dt. (21) 

J 0 Jo 


Proof: We only give the proof of the second equa¬ 
tion since the proof of the first one is identical. From 
the assumptions in the statement of the lemma, we have that 
z(t) > 0 for every t > 0 with probability one. Therefore, by 
the linearity of expectations, the identity ||v||i = l„v that 
holds for a general v £ R7, and the fact that = 1, we 
can show that t T Mr E[^(t) <g> z (t)} = E[(t r M £,(t))(t T r z{t))] = 
ljj£[z(i)] = ||£7[z(£)]|| 1 . Integrating the both sides of this 
equation with respect to t from 0 to 00 , we obtain GD- ■ 


2 A preliminary version of the theorem can be found in |14) . Also, we 
remark that the theorem is a continuous-time counterpart of 1251 Theorem 2]. 


We now have the elements needed to prove Theorem [ 4 ] 
Proof of Theorem [7} First assume that E is internally 
mean stable and ||E|| j < 7. Then, A is Hurwitz stable by 
Proposition IT] and, therefore, invertible. Thus, the vector 77 = 
1 ~Mr(D — C'A~ 1 B) — 7Iis well-defined. Let us first 
show ?7 < 0 . Take e > 0 such that ||£||i < 7 — e. Then, 
for every initial state op and w £ £i(R + ,R+), if xo = 0 , 
then E[z] £ £i(R+,R© and ||£[z]|| £l < ^7 - e)||w|| £l . 
Therefore, by Lemma [3J we have that l~l Ar . J^° E[£ ® z ]dt < 
(7 — e)lj 4s /“(£[£] ® w) dt. Then, by ( fl 9 | ), we obtain 

poo 

[r] + el m s ) / (E[£\ ® w )dt < 0. (22) 

Jo 

Now, let i £ [M] and j £ [s] be arbitrary. Let ej and fj be the 
*-th and j-th standard unit vectors in R M and R s , respectively. 
Let (Jo = i and, for every r > 0, define w T = rx[o,i/r]fj £ 
£i(R+,R+), where X[o,i/r] : f {0,1} denotes the indica¬ 

tor function of the interval [0,1/r]. Notice that, by a stan¬ 
dard argument in distribution theory |24j, the function r i—>■ 
t X[o,i/t] converges to the Dirac delta function as r —> oo in 
the space of distributions. Thus, in the limit of r —► oo, we ob¬ 
tain J n °°(£[£] <g> w r ) dt —j (£[£(0)]) ® fj = ej g> fj. There¬ 
fore, |22] shows that (77 + el~M S )( e i ® fj) £ 0- Since % and j 
are arbitrary, it must be the case that i)+el]^ s < 0 and, hence, 
r/ < 0, as we wanted to show. Now, by Proposition [2] there 
exists a positive v £ R Mn such that 

v T A + lj[)f r C <0, v T B -(- ±M r D < (23) 

Then, from the definition of the matrix A (see Proposition^, it 
is easy to see that the inequalities ( fT7| i are satisfied by the pos¬ 
itive vectors vi,..., vju £ R” given by v = (vj, ..., 

This completes the proof of the necessity part of the theorem. 

On the other hand, assume that there exist positive vec¬ 
tors Vi,...,vm £ R" such that ( |TT| holds. Then we can 
see that the positive vector v = (v^, ■ • •, vJf) T £ R M " 
satisfies Therefore, by Proposition [2] the linear time- 

invariant positive system (A, B,C, D) is stable and has the 
£i-gain less than 7. Thus, A is Hurwitz stable and, by 
Proposition [T] the Markov jump linear system E is internally 
mean stable. We need to show that ||£||i < 7. Let op and 
w £ £i(R + , R^_) be arbitrary and set xo = 0. Let z denote the 
corresponding trajectory of E. From ( |2T| and (p~9] ), it follows 
that / 0 °°||£[z(f)]|| 1 dt = lJf r (D — CA~ 1 B) / 0 (£[(]®w) dt. 
Since the positive linear time-invariant system (A, B,C, D) 
has the £i-gain less than 7, Proposition [2] shows that there 
exists an e > 0 satisfying 1 T Mr {D — CA~^B) < (7 — e)l r Ms . 
Therefore, we obtain / 0 3 °||£[z(f)]||i dt < (7— ejjlwHpj, where 
we used ( ]20| . This proves E[ z] £ £i(R + ,R^_) and therefore 
the Ex-stability of E. Moreover we have ||E||i < 7 — e < 7. 
This completes the proof of the theorem. ■ 

Now we can give the proof of Theorem [3] 

Proof of Theorem jlj From Theorem [4]the Performance- 
Constrained Disturbance Attenuation problem can be stated as 
the following optimization problem: 

minimize R 

^,{Fk,Gk,Hk}k£V A v i}i£[M] 

subject to Vi > 0, and <0- 



Applying the change of variables in m we obtain the 
optimization program stated in the theorem. We can also see 
that the optimization problem is, in fact, a matrix geometric 
program by following the explanation in Remark |3| ■ 

Remark 6: As mentioned in Subsection |III-C| Problem II-B 
considers the case of minimizing 7 (i.e., maximizing the dis¬ 
turbance attenuation) under a budget constraint. The statement 
and solution of Problem II-B fall straightforward from the 
cases previously considered and, hence, details are omitted. 

VI. Numerical Simulations 
In this section, we illustrate our network design frame¬ 
work to stabilize the dynamics of a disease spreading in a 
time-varying network of individuals. We consider a popular 
networked dynamic model from the epidemiological liter¬ 
ature, the networked Susceptible-Infected-Susceptible (SIS) 
model |23[ . According to this model, the evolution of the 
disease in a networked population can be described as: 

ik(t) = -6 k x k (t ) + I3 k ^2 K ke (t)xe(t) + e k w k {t) (24) 
aev 

for k £ V, where x k (t) is a scalar variable representing the 
probability that node k is infected at time t. The parame¬ 
ter S k > 0, called the recovery rate , indicates the rate at 
which node k would be cured from a potential infection. The 
parameter /3 k > 0, called the infection rate, indicates the rate 
at which the infection is transmitted to node k from its infected 
neighbors. The exogenous signal e k w k (t), where e k > 0 is a 
constant and w k is an R + -valued function, is introduced to 
explain possible transmission of infection from outside of the 
network. The entries of the time-varying adjacency matrix of 
the contact network are K k ^{t) £ {0,1}, for k,£ £ V. 

We consider the following epidemiological problem & 
Assume we have access to vaccines that can be used to reduce 
the infection rates of individuals in the network, as well as 
antidotes that can be used to increase their recovery rates. 
Assuming that both vaccines and antidotes have an associ¬ 
ated cost, how would you distribute vaccines and antidotes 
throughout the individuals in the network in order to eradicate 
an epidemic outbreak at a given exponential decay rate while 
minimizing the total cost? We state this question in rigorous 
terms below and present an optimal solution using geometric 
programming. Let C\(f3 k ) and 02 ( 8 k ) denote the costs of tuning 
the infection rate /3 k and the recovery rate 5 k of agent k, 
respectively. We assume that these rates can be tuned within 
the following feasible intervals: 

0 < @ k < fi k < j3 k , 0<S k <5 k < S k ■ (25) 

The problem of finding the optimal allocation of vaccines 
and antidotes in a static network was recently solved in [ r9j|, 
under certain assumptions on the cost functions ci and c 2 . In 
what follows, we solve this problem for time-varying contact 
networks, K(t) = K a r t \, where 0 is a time-homogeneous 
Markov process. In particular, when there is no disturbance 
from outside of the network, i.e., when e k = 0 for every 
k £ V, we can formulate this problem as a particular version 
of the Performance-Constrained Stabilization problem solved 
in Section [IV] as follows: 


Problem 1: Assume e k = 0 for every k £ V. Given a 
desired decay rate A > 0, tune the spreading and recovery 
rates {/ 3 k } k& y and {<5fc}fc e v in the network such that the 
disease modeled in ( [24] ) is eradicated at an exponential decay 
rate of A and a minimum cost R = Efcev( Cl (^ fc ) + c 2 {Sk)), 
while satisfying the box constraints in ( |25j ). 

We can transform the above problem into a Performance- 
Constrained Stabilization problem, as follows. It is easy to 
see that the systems ( [24] ) form a network of positive linear 
systems x k = —5 k x k + /3 k u k + e k w k with the coupling 
Uk(t) = Er e v K k n(t)p k x t (t). We set f k (F k ,G k ,H k ) = 
C 2 (—F k ) for each k £ V and <f>( k ,e) = c i/d k for each 
(k,£) £ £, where d k is the indegree of vertex k defined by 
d k = \{(. £ V : (k,£) £ £}\. For illustration purposes, we use 


the set of cost functions proposed in 119]: 


ci(/3fc) = 


ft 1 -?: 1 




—, c 2 (<5fc) = 


{l - Sk )- 1 - (l - Sk )- 1 

(1 - ^fc)- 1 - (1 - 1J*;)- 1 ' 


Notice that c\ is decreasing, c 2 is increasing, and the range 
of ci and c 2 are both [0,1], In this particular example, we let 
A k = 1 for every k £ V. Then, it follows that 


f k (F k ,G k , H k ) 


(1 - Jfc )- 1 - (1 - ' 


The negative constant term — (1 — &) -1 in the numerator of 
f k can be ignored, since it only changes the value of the total 
cost R by a fixed constant (while the optimal allocation is 
unchanged). For the same reason, we can ignore the negative 
constant — (3 k in the numerator of (j>r k ,e)- Notice that, after 
ignoring these constant terms, the functions f k and 4 >{k,t) are 
posynomials and, hence. Assumption [3] holds true. Notice also 
that, since fj k does not depend on l, we need the additional 
constraints = ■ ■ • = T( kt m, which can be implemented 

as monomial constraints T( fc = 1 for f = 2 ,N. 

Once the cost functions are defined, let us consider a 
particular model of network switching to illustrate our results. 
Consider a set of agents 1 ,,N divided into h disjoint 
subsets called Households. Each Household has exactly one 
agent called a Worker. The set of Workers is further divided 
into w disjoint subsets called Workplaces. We assume that 
the topology of the contact network switches between three 
possible graphs. In the first contact graph K\, a pair of agents 
are adjacent if they are in the same household. The second 
graph K '2 is a random graph in which any pair of Workers 
are adjacent with probability p = 0.3, independently of other 
interactions. We also assume that non-Worker agents are adja¬ 
cent if they belong to the same Household. The third graph A 3 
is a collection of disconnected h + w complete subgraphs. In 
particular, we have w complete subgraphs formed by disjoint 
sets of Workers sharing the same Workplace, and h complete 
subgraphs formed by sets of non-Workers belonging to the 
same Household. 

We consider the case in which the topology of the network 
switches according to a Markov process cr taking its values in 
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Fig. 1: Correction versus prevention per node. Circle: Workers. 
Square: non-Workers. The color of markers indicate the size of 
Workplaces/Households where agents belong to. 


the set {1,2, 3, 2} with the following infinitesimal generator: 
'-1/13 1/13 0 0 


n = 


-l 

o 

o 


i 

-1/9 

0 


0 

1/9 

-1 


(26) 


This choice reflects Workers’ simplified schedule: 13 hours of 
stay at home (6 PM-7 AM), 1 hours of commute (7 AM-8 
AM), 9 hours of work (8 AM-5 PM), and 1 hour of commute 
(5 PM-6 PM). 

In this setup, we find the optimal allocation of resources 
to control an epidemic outbreak for the following parameters. 
First, we randomly generate a set of n = 247 agents with h = 
71 Workers (as many as Households) and w = 10 Workplaces. 
We let 4 = 0.01, 4 = 0.05, 4 = -0.5, 4 = -0.1, and 
efc = 0 for every k € V. We then solve the Performance- 
Constrained Stabilization problem (Problem [I]) with A = 0.01 
using Theorem |T| and obtain the optimal investment strategy 
with total cost of R = 54.38. 

Fig. □ is a scatter plot showing the relationship between 
the investments on 4 (corrective actions) and 4 (preventive 
actions) for each agent. We can see that, in order to maximize 
the effectiveness of our budget, we need to invest heavily 
on Workers, in particular on those belonging to larger Work¬ 


places. Figs. 2a and 2b show the infection probabilities Xk 
starting from t = to = 8, when Workers start working 
(therefore a (to) = 3). We choose the vector of initial proba¬ 
bilities of infection, Xk(to), at random from the set {0,1} 4 
Each solid curve in the figures represents the evolution in the 
probability of infection of an agent over time. The dashed 
vertical lines indicate the times when the graph changes. 
We are also including a dashed curve representing £ : [||x(f)||] 
when a follows the Markov process with the infinitesimal 
generator |26|. 

In the second example, we randomly set ej, to be ei¬ 
ther 0 or 1 for each k £ V, and solve the Performance- 
Constrained Disturbance Attenuation problem with the con¬ 
straint sup^ ez:i(R+;M « ) (|| J E , [x]|| z:i /||w|| z:i ) < 7 = 40. We 
use the same cost functions as in our first example. Using 
Theorem [3] we obtain the optimal values of 4 and 4 with 
a total cost of 67.81. Fig. [3] is a scatter plot showing the 



(a) Infection probabilities of Workers. 



(b) Infection probabilities of non-Workers. 

Fig. 2: Solid: Infection probabilities Xk(t) of Workers |a| and non- 
Workers (Jb]». Dotted: 25[||a:(t)||i]. The colors of the lines indicate the 
size of Workplaces/Households where agents belong to. The labels 
’W’ and "H’ indicate the first two time-periods when Workers are 
in their Workplaces and Households, respectively. Narrow vertical 
stripes represent the time-periods when Workers are commuting. 


relationship between the investments on 4 and 4- We can 
observe different patterns of resource allocation for agents with 
and without disturbance; in general, those with disturbance 
receive more allocation for corrective resources, while those 
without disturbance do more for preventive resources. 

Before closing this section, we briefly discuss the computa¬ 
tional cost of solving the geometric programs in this example. 
For simplicity, we assume that the number of nonzero entries 
in each column of n is less than or equal to N. Then, 
according to the notation in Proposition [3] we can show that 
n = N(M + 2), k < 2 N, and m < N(M + 4) for Problem [I] 
Therefore, by Proposition [3] the computational cost for solving 
Problem [l] equals 0(N 7 /^M 7 / 2 ). Similarly, we can show that 
the same computational cost is required by the Performance- 
Constrained Disturbance Attenuation problem. 

VII. Conclusions and Discussion 

In this paper, we have proposed an optimization framework 
to design the subsystems and coupling elements of a time- 
varying network to satisfy certain structural and functional 
requirements. We have assumed there are both implementation 
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Fig. 3: Correction versus prevention per node. Circle: Workers. 
Diamond: Workers with disturbance. Square: non-Workers. Triangle: 
non-Workers with disturbance. The color of markers indicate the size 
of Workplaces/Households where agents belong to. 


costs and feasibility constraints associated with these network 
elements, which we have modeled using posynomial cost 
functions and inequalities, respectively. In this context, we 
have studied several design problems aiming at finding the 
cost-optimal network design satisfying certain budget and 
performance constraints, in particular, stabilization rate and 
disturbance attenuation. We have developed new theoretical 
tools to cast these design problems into geometric programs 
and illustrated our approach by solving the problem of sta¬ 
bilizing a viral spreading process in a time-switching contact 
network. 

A possible direction for future research is the extension to 
networks of nonpositive systems. The proposed framework 
is not directly applicable to this case due to the positivity 
constraint in geometric program. However, by applying the 
framework to the upper-bounding linear dynamics of vector 
Lyapunov functions presented in ED’ we may be still able 
to accomplish an efficient, if not optimal, network design for 
stabilization and disturbance attenuation. 
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